In [1]:
%matplotlib inline
from matplotlib import pyplot as plt
from os import path
import numpy as np
import cv2
import pandas as pd
def show_images(images,titles=None, scale=1.3):
"""Display a list of images"""
n_ims = len(images)
if titles is None: titles = ['(%d)' % i for i in range(1,n_ims + 1)]
fig = plt.figure()
n = 1
for image,title in zip(images,titles):
a = fig.add_subplot(1,n_ims,n) # Make subplot
if image.ndim == 2: # Is image grayscale?
plt.imshow(image)
else:
plt.imshow(cv2.cvtColor(image, cv2.COLOR_RGB2BGR))
a.set_title(title)
plt.axis("off")
n += 1
fig.set_size_inches(np.array(fig.get_size_inches(), dtype=np.float) * n_ims / scale)
plt.show()
The images supplied vary drastically in their palettes. To simplify feature extraction we need to first normalize them to the same color.
In [2]:
img_path = "/kaggle/retina/train/sample"
# this is the image into the colors of which we want to map
ref_image = cv2.imread(path.join(img_path, "6535_left.jpeg"))
# images picked to illustrate different problems arising during algorithm application
image_names = ["16_left.jpeg", "10130_right.jpeg", "21118_left.jpeg"]
image_paths = map(lambda t: path.join(img_path, t), image_names)
images = map(lambda p: cv2.imread(p), image_paths)
# let's see what we've got
image_titles = map(lambda i: path.splitext(i)[0], image_names)
show_images(images, image_titles, scale = 0.8)
show_images([ref_image], ["ref"], scale = 0.8)
In [3]:
#Pyramid Down & blurr
# Easy-peesy
def pyr_blurr(image):
return cv2.GaussianBlur(cv2.pyrDown(image), (7, 7), 30.)
ref_image = pyr_blurr(ref_image)
images = map(lambda i: pyr_blurr(i), images)
In [4]:
def display_contours(image, contours, color = (255, 0, 0), thickness = -1, title = None):
imShow = image.copy()
for i in range(0, len(contours)):
cv2.drawContours(imShow, contours, i, color, thickness)
show_images([imShow], scale=0.7, titles=title)
image = ref_image
# this is a good threshold for Canny edge finder, but it does not always work. We will see how to deal with it furhter on.
thresh = 4
# Searcing for the eye
# Let's see how this works setp-by-step
# convert to a one channel image
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# Canny edge finder
edges = np.array([])
edges = cv2.Canny(gray, thresh, thresh * 3, edges)
# Find contours
# second output is hierarchy - we are not interested in it.
contours, _ = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
# Let's see what we've got:
display_contours(image, contours, thickness=2)
print "{:d} points".format(len(np.vstack(np.array(contours))))
In [5]:
# Now let's get only what we need out of it
hull_contours = cv2.convexHull(np.vstack(np.array(contours)))
hull = np.vstack(hull_contours)
# we only get one contour out of it, let's see it
display_contours(image, [hull], thickness=2, color=(0, 255, 0))
In [6]:
# Now let's create a mask for this image
def createMask((rows, cols), hull):
# black image
mask = np.zeros((rows, cols), dtype=np.uint8)
# blit our contours onto it in white color
cv2.drawContours(mask, [hull], 0, 255, -1)
return mask
mask = createMask(image.shape[0:2], hull)
show_images([mask])
In [7]:
def find_eye(image, thresh = 4, size=256):
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# Canny edge finder
edges = np.array([])
edges = cv2.Canny(gray, thresh, thresh * 3, edges)
# Find contours
# second output is hierarchy - we are not interested in it.
contours, _ = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
# Now let's get only what we need out of it
hull_contours = cv2.convexHull(np.vstack(np.array(contours)))
hull = np.vstack(hull_contours)
def createMask((rows, cols), hull):
# black image
mask = np.zeros((rows, cols), dtype=np.uint8)
# blit our contours onto it in white color
cv2.drawContours(mask, [hull], 0, 255, -1)
return mask
mask = createMask(image.shape[0:2], hull)
# returning the hull to illustrate a few issues below
return mask, hull
In [8]:
# here is the histogram of the reference image:
hist_green = cv2.calcHist([ref_image], [1], None, [256], np.array([0, 255]))
fig = plt.figure()
plt.bar(np.arange(256), hist_green, width=2, color='g', edgecolor='none')
fig.set_size_inches(np.array(fig.get_size_inches(), dtype=np.float) * 2)
plt.show()
In [9]:
# we should really take it over the masked region, to remove the irrelevent background
mask, hull = find_eye(ref_image)
hist_green = cv2.calcHist([ref_image], [1], mask, [256], np.array([0, 255]))
fig = plt.figure()
plt.bar(np.arange(256), hist_green, width=2, color='g', edgecolor='none')
fig.set_size_inches(np.array(fig.get_size_inches(), dtype=np.float) * 2)
plt.show()
In [10]:
old_images = images
images = [ref_image]
def saturate (v):
return map(lambda a: min(max(round(a), 0), 255), v)
def get_masks(images):
return map(lambda i: find_eye(i)[0], images)
maskHull = get_masks(images)
channels = map(lambda i: cv2.split(i), images)
imMask = zip(channels, maskHull)
nonZeros = map(lambda m: cv2.countNonZero(m), maskHull)
# grab three histograms - one for each channel
histPerChannel = map(lambda (c, mask): \
[cv2.calcHist([cimage], [0], mask, [256], np.array([0, 255])) for cimage in c], imMask)
# compute the cdf's.
# they are normalized & saturated: values over 255 are cut off.
cdfPerChannel = map(lambda (hChan, nz): \
[saturate(np.cumsum(h) * 255.0 / nz) for h in hChan], zip(histPerChannel, nonZeros))
# let's see the cdf over the green channel
fig = plt.figure()
plt.bar(np.arange(256), cdfPerChannel[0][1], width=2, color='g', edgecolor='none')
fig.set_size_inches(np.array(fig.get_size_inches(), dtype=np.float) * 2)
plt.show()
In [11]:
# now compute the map. We will need the reference image cdf for that, so we might as well put it all together:
def saturate (v):
return np.array(map(lambda a: min(max(round(a), 0), 255), v))
def get_masks(images):
return map(lambda i: find_eye(i)[0], images)
def calc_hist(images, masks):
channels = map(lambda i: cv2.split(i), images)
imMask = zip(channels, masks)
nonZeros = map(lambda m: cv2.countNonZero(m), masks)
# grab three histograms - one for each channel
histPerChannel = map(lambda (c, mask): \
[cv2.calcHist([cimage], [0], mask, [256], np.array([0, 255])) for cimage in c], imMask)
# compute the cdf's.
# they are normalized & saturated: values over 255 are cut off.
cdfPerChannel = map(lambda (hChan, nz): \
[saturate(np.cumsum(h) * 255.0 / nz) for h in hChan], \
zip(histPerChannel, nonZeros))
return np.array(cdfPerChannel)
# compute color map based on minimal distances beteen cdf values of ref and input images
def getMin (ref, img):
l = [np.argmin(np.abs(ref - i)) for i in img]
return np.array(l)
# compute and apply color map on all channels of the image
def map_image(image, refHist, imageHist):
# each of the arguments contains histograms over 3 channels
mp = np.array([getMin(r, i) for (r, i) in zip(refHist, imageHist)])
channels = np.array(cv2.split(image))
mappedChannels = np.array([mp[i,channels[i]] for i in range(0, 3)])
return cv2.merge(mappedChannels).astype(np.uint8)
# compute the histograms on all three channels for all images
def histogram_specification(ref, images, masks):
cdfs = calc_hist(images, masks)
mapped = [map_image(images[i], ref[0], cdfs[i, :, :]) for i in range(len(images))]
return mapped
In [12]:
# the fruit of our labor
images = old_images
refMask = get_masks([ref_image])
refHist = calc_hist([ref_image], refMask)
masks = get_masks(images)
histSpec = histogram_specification(refHist, images, masks)
print "Reference Image"
show_images([ref_image], ["Ref"], scale = 0.9)
print "Original Images"
show_images(images, titles = image_titles, scale = 0.9)
print "Mapped Images"
show_images(histSpec, titles = image_titles, scale = 0.9)
At this point it seems like a good idea to filter out the background noise which appears to be present in image 3: to at least make sure the background is entirely black, so that future steps of the algorithm could mask it out easily.
Not a problem:
In [13]:
def mask_background(image, mask):
channels = np.array(cv2.split(image))
# now mask off the backgrownd
for i in range(0, 3):
channels[i, mask == 0] = 0
return cv2.merge(channels).astype(np.uint8)
maskedBg = [mask_background(h, m) for (h, m) in zip(histSpec, masks)]
print "Masked Background"
show_images(maskedBg, titles=image_titles, scale=0.9)
Now what?! This does not look desirable at all! Since we were histogramming using a mask, and our mask appears to have consumed so much of the actual image, we need to do something about it. Let's take a look at the mask:
In [14]:
show_images([masks[2]], scale = 0.9)
Right. This is not the mask we are looking for. As we can see, the image is too dark, so to retain more information we need to lower the threshold. A simple heuristic works pretty well: if the resulting mask occupies less than 41% percent of the image, we reduce the threshold to 1. We must be careful, though. If we approach the problem with this threshold at the very beginning i.e, set the threshold to 1 for all images, for many images we will catch a lot of background noise in our mask and it will be no mask at all! (the reverse effect).
In [22]:
mask, _ = find_eye(images[2], 1)
show_images([mask], scale = 0.9)
let's apply that:
In [15]:
mask, _ = find_eye(images[2], 4)
masks[2] = mask
histSpec = histogram_specification(refHist, [images[2]], [masks[2]])
print "Defective mask:"
show_images(histSpec, titles = image_titles[2:], scale = 0.9)
mask, _ = find_eye(images[2], 1)
masks[2] = mask
histSpec = histogram_specification(refHist, [images[2]], [masks[2]])
print "Corrected Mask:"
show_images(histSpec, titles = image_titles[2:], scale = 0.9)
maskedBg = [mask_background(h, m) for (h, m) in zip(histSpec, [masks[2]])]
print "Masked Background"
show_images(maskedBg, titles=image_titles[2:], scale=0.9)
And a brief illustration of what happens if we set the threshold too low from the start:
In [16]:
#Mask that catches too much:
mask, _ = find_eye(images[1], 1)
show_images([mask])
masks[1] = mask
histSpec = histogram_specification(refHist, [images[1]], [masks[1]])
print "Defective mask:"
show_images(histSpec, titles = image_titles[1:2], scale = 0.9)
# Reset the threshold to 4 and re-run
mask, _ = find_eye(images[1], 4)
masks[1] = mask
histSpec = histogram_specification(refHist, [images[1]], [masks[1]])
print "Corrected Mask:"
show_images(histSpec, titles = image_titles[1:], scale = 0.9)
maskedBg = [mask_background(h, m) for (h, m) in zip(histSpec, [masks[1]])]
print "Masked Background"
show_images(maskedBg, titles=image_titles[1:], scale=0.9)